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1.0  INTRODUCTION 


( 

The  Performance  Modeling  of  Autonomous  Electro-Optical  Sensors  (PM)  program  ia  a 
basic  research  effort  aimed  at  the  fundamental  analysis  of  image  processing  methodology. 
To  this  end,  an  investigation  into  the  the  use  of  LAD  AR  imagery  for  target  recognition  was 
initiated.  Preliminary  research  and  methodology  development  occurred  from  November 
1986  to  May  1987.  Research  and  implementation  are  continuing  and  this  report  details 
the  investigation  performed  and  results  obtained  during  the  period  from  May  15,  1987 
through  November  15,  1987. 

During  the  previous  six  months,  the  LADAR  target  recognition  investigation  has  fo¬ 
cused  on  the  area  of  object  segmentation.  The  issues  that  were  addressed  include: 

Si  Surface  segmentation  via  region  growing^ 

4  Robust  planar  surface  estimation , 

Sr  Background  plane  removal  , 

S  •  Smoothing  of  the  image  prior  to  the  estimation  of  derivatives  ^  #  J 
'"--m  Crease  edge  detection  / 

The  introduction  to  Section  2.0  contains  a  general  problem  definition  and  a  brief  method¬ 
ology  overview.  The  various  subsections  describe  those  aspects  of  segmentation  that  were 
addressed  in  this  reporting  period.  Specific  techniques  being  investigated  are  detailed, 
along  with  any  pertinent  demonstrations  that  have  been  performed.  More  information 
on  the  overall  recognition  methodology  and  matching  in  particular  afr^centained  in  the 
Second  Semi-Annual  Report  (7).  Areas  of  future  work  and  investigation  ard^mmmarized 
in  Section  3.0. 


2.0  TECHNICAL  PROGRESS 

la  tb«  Second  Semi-Annual  Report  [7]  we  described  our  approach  to  target  recognition 
based  upon  laser  radar  sensor  information.  This  approach  emphasizes  (i)  the  3-dimensional 
geometry  of  the  targets  and  the  scene,  (ii)  both  local  and  global  features,  and  (iii)  model- 
based  graph  matching.  In  [7]  an  overview  of  segmentation,  description,  and  matching 
were  given,  and  detection  was  touched  upon  briefly.  The  major  focus  of  investigation 
during  that  period  involved  matching,  with  a  proof  of  concept  example  demonstrating  the 
feasibility  of  the  selected  approach.  During  the  last  six  months,  the  focus  has  shifted  to 
the  development  of  an  effective  segmentation  methodology. 

This  report  contains  a  detailed  discussion  of  the  areas  that  were  investigated  during  the 
segmentation  study.  Overviews  of  the  algorithms  that  are  being  researched  are  presented 
along  with  examples  of  the  associated  experimental  results. 

2.1  Segmentation/Description  Overview 

With  range  imagery  the  segmentation  step  involves  determining  the  extended  structure 
of  an  object  within  the  given  target  window.  Specifically,  this  procedure  entails:  (i)  the 
detection  of  target  boundaries,  (ii)  the  location  of  edges  between  surfaces  of  the  target, 
and  (iii)  the  determination  of  target  surfaces,  but  not  necessarily  in  that  order. 

Segmentation  of  intensity  images  is  a  well  studied  problem  in  image  understanding 
and  computer  vision  applications.  On  the  contrary,  segmentation  of  range  images  for 
target/ background  separation,  and  for  detecting  edges  between  surfaces  is  a  relatively  new 
research  topic. 

There  are  two  distinct  classes  of  methods  for  detecting  edges.  The  first  is  to  search  for 
edges  directly  by  hunting  for  properties  that  distinguish  an  edge.  There  axe  many  different 
edge  detectors  of  this  type,  some  of  which  are  more  robust  to  noise  than  others.  Another 
method  for  finding  edges  is  indirect.  First,  the  surfaces  of  the  object  are  identified  and 
then  the  intersection  of  two  adjacent  surfaces  indicates  the  presence  of  an  edge.  Through 
the  combination  of  these  two  approaches,  more  accurate  estimates  of  both  surface  and 
edge  information  can  be  obtained. 

An  iterative  region  growing  technique  for  identifying  surfaces  also  provides  a  “para 
metric”  representation  of  the  surfaces .  Such  a  represent,  tion  is  critical  for  effective  object 
description.  Specifically,  surface  parameterization  is  reqi  red  to  reduce  the  description  of 
a  geometric  entity  to  a  small  set  of  measures  which  are  invariant  to  rotations,  translations, 
and  changes  in  scale.  In  other  words,  it  is  he  first  step  in  forming  a  general  object  descrip¬ 
tion  that  can  be  employed  for  matching.  M  ’  -re  by  sing  a  surface  based  segmentation 
in  conjunction  with  edge  detection  the  fun  lion,  of  segmentation  and  initial  description 
can  be  accomplished  simultaneously  with  a  high  degree  of  accuracy. 

Besl  and  Jain  jl|  have  developed  a  fast,  region  growing  segmentation  algorithm  which 
starts  from  seed  regions  based  on  uniform  invariant  curvature  classification  and  proceed-, 
hy  parallel,  iterative  region  growing.  We  are  using  this  surface  based  algorithm  as  the 
core  of  our  segmentation  approach.  Section  2.2  describes  this  technique  along  with  our 
modifications,  and  presents  .some  examples  of  its  behavior. 
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One  of  the  most  important  components  of  the  region  growing  algorithm  is  the  surface 
estimation  which  is  effectively  reduced  to  planar  estimation.  During  preliminary  tests  of 
the  matching  technique  in  February  1987  [4],  it  was  observed  that  planar  estimation  of 
the  noisy  range  data  was  extremely  erratic  when  a  least  square  fit  was  used.  Therefore 
a  comparison  study  of  various  planar  estimation  techniques  was  performed.  Section  2,3 
describes  the  various  estimation  methods  that  were  compared  as  well  as  the  tests  that  were 
performed  and  their  results. 

Another  issue  which  arose  was  the  appropriate  estimation  of  thresholds  used  during 
the  region  growing  process.  Many  of  the  thresholds  are  computed  from  the  image,  but 
may  be  set  incorrectly  if  too  much  background  is  included  in  this  computation.  This 
occurs  because  the  background  pixels  often  have  a  larger  variance  from  a  planar  tit  since 
(i)  they  generally  don’t  come  from  exact  planes  and  (ii)  a  high  angle  of  incidence  makes 
the  footprint  of  the  pixel  cover  a  wide  interval  of  ranges.  Ideally  the  thresholds  should  be 
calculated  only  from  pixels  on  the  target  itself.  Using  LADAR  imagery  alone  though  it 
is  not  easy  to  extract  the  target  boundary.  Section  2.4  discusses  the  problems  associated 
with  finding  and  eliminating  a  large  percentage  of  background  pixels  during  pre-processing 
for  the  region  grower. 

During  the  region  growing,  and  for  most  types  of  edge  detection,  it  is  necessary  to  esti¬ 
mate  various  partial  derivatives  from  the  image.  Due  to  noise  and  quantization  error  in  the 
image,  these  estimates  are  often  unreliable.  Smoothing  can  be  used  to  improve  derivative 
estimation  if  it  is  performed  appropriately.  Section  2.5  discusses  how  regularization  can 
be  employed  to  optimally  filter  an  image  for  a  specified  operation. 

The  two  types  of  edges  that  occur  in  range  imagery  are  (i)  jump  edges  which  are  formed 
when  there  is  a  discontinuity  (other  than  a  range  ambiguity)  in  the  relative  range,  and  (ii) 
crease  edges,  which  are  derived  from  discontinuities  in  orientation  without  discontinuities  in 
range.  Jump  edges  can  be  identified  quite  accurately  either  directly,  or  indirectly  using  the 
region  growing  algorithm.  Crease  edges  can  be  more  difficult  to  identify  via  region  growing 
and  their  location  may  vary.  Therefore  an  alternate  crease  edge  detection  technique  has 
been  implemented  to  supplement  and  correct  the  output  of  the  surface-based  segmentation. 
The  crease  edge  detection  technique  used  is  a  variant  of  that  proposed  by  Mitiche  and 
Aggarwal  in  |8j.  This  technique  is  described  in  section  2.6  and  an  example  is  presented 
demonstrating  its  function. 

2.2  Surface-Based  Segmentation 

The  goal  of  surface-based  segmentation  is  to  accurately  estimate  the  true  underlying  surface 
shape  of  objects  in  LADAR  range  imagery.  The  technique  developed  by  Bes)  and  Jain  [ij 
employs  iterative  region  growing,  and  is  based  upon  two  significant  assumptions  regarding 
the  characteristics  of  targets  in  range  data. 

*  All  objects  of  interest  may  be  represented  by  piecewise  smooth  surfaces. 

•  Planar  and/or  quadratic  surfaces  are  sufficient  to  mode!  targets  under  typical  con¬ 
ditions  of  sensor  resolution. 


Seed  regions  form  the  starting  point  for  the  region  growing  process.  They  should 
represent  small,  homogeneous  regions  not  containing  any  boundary  points  and  be  a  reliable 
representative  of  the  surrounding  surface  curvature.  Beni  and  Jain  use  the  Gaussian  ( K ) 
and  mean  (H)  curvature  to  form  a  HK-Sign  map  which  divides  the  image  into  regions 
belonging  to  an  invariant  curvature  class.  The  eight  possible  curvature  classes  -  peak,  pit, 
ridge,  valley,  flat,  minimal,  saddle  valley,  and  saddle  ridge  -  are  used  as  a  starting  point 
when  forming  seed  regions.  Gaussian  and  mean  curvature  are  defined  and  discussed  in 
more  detail  in  [1,7). 

A  summary  of  the  major  steps  involved  in  this  algorithm  is  presented  here. 

1.  Divide  the  image  into  regions  of  invariant  surface  class.  This  is  done  using  the 
HK-Sign  map. 

2.  Take  the  largest  connected  region  of  invariant  surface  class  and  contract  it  so  that 
only  points  interior  to  the  region  remain.  This  is  the  seed  region. 

3.  Fit  a  plane  (or  higher  order  surface)  to  the  seed  region.  This  is  the  estimated  surface. 

4.  Search  for  all  the  points  in  the  image  that  are  approximately  fit  by  the  estimated 
surface  and  find  the  largest  connected  subset  which  overlaps  the  seed  region.  This 
becomes  the  new  seed  region. 

5.  Check  conditions  for  stopping.  If  region  growing  is  to  continue  go  bacl  to  step  3. 
There  are  a  variety  of  conditions  related  to  stopping  at  this  point.  Briefly,  the  more 
important  checks  involve  (i)  the  error  in  the  surface  fit  vs.  a  calculated  threshold,  (ii) 
the  relative  change  in  region  size  since  the  previous  iteration,  and  (iii)  the  current 
number  of  iterations. 

6.  If  the  final  region  is  acceptable  then  add  it  to  the  surface  list  and  remove  the  points 
in  that  region  from  the  surface  class  image.  Otherwise,  remove  the  original  seed 
region  points  from  the  surface  class  image.  This  is  done  to  prevent  looping  caused 
by  repeatedly  trying  to  fit  the  same  surfaces.  The  acceptance  of  a  region  is  also  based 
upon  a  variety  of  conditions,  related  largely  to  region  size,  error  of  fit,  and  whether 
this  region  is  contained  within  a  previously  accepted  region, 

7.  Go  back  to  choose  another  surface  class  region  in  step  2. 

The  purpose  of  this  investigation  was  to  determine  whether  this  technique  could  be 
effectively  used  for  segmenting  surfaces  of  vehicular  targets.  To  address  this  question 
several  alterations  were  made  to  Besl’s  algorithm. 

*  initial  runs  with  the  algorithm  showed  a  tendency  to  use  a  parabolic  fit.  over  the 
entire  target.  Since  the  target  types  of  interest  are  composed  of  planar  surfaces,  the 
surface  fitting  was  restricted  to  planes.  This  restriction  is  being  used  for  these  initial 
tests  and  may  be  relaxed  later. 
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®  The  planar  surface  fitting  is  being  performed  using  an  estimation  technique  that  our 
analysis  showed  to  be  superior  to  least  squares  under  contaminated  noisy  conditions. 
This  technique  is  called  M-estuxtafcefl.  It  is  presented  in  section  2.3  along  with  the 
tests  comparing  it  to  alternative  estimation  methods. 

•  With  the  addition  of  M-estim&tes  for  plane  fitting,  a  portion  of  the  code  which 
attempted  to  remove  outliers  in  a  very  adhoc  manner  was  removed.  M-estim&tes 
down-weights  the  real  outliers  much  more  precisely. 

•  The  variance  of  the  gradient  in  the  background  was  often  higher  than  on  the  target 
itself,  and  many  of  the  thresholds  are  set  based  upon  the  perceived  gradient  variance 
in  the  entire  image.  Therefore,  the  thresholds  were  being  skewed  due  to  background 
noise  and  the  segmentation  of  subregions  within  a  target  was  poor  or  non-existent.  By 
masking  out  a  sizable  fraction  of  the  background  pixels  when  calculating  thresholds 
the  performance  was  much  improved.  Since  surface  segmentation  generally  occurs 
after  the  target  is  mostly  separated  from  the  background,  a  mask  is  very  appropriate. 
Section  2.4  discusses  how  the  background  masks  can  be  obtained. 

•  The  region  growing  algorithm  does  not  take  into  consideration  any  explicit  edge 
information.  Jump  edges  are  usually  not  a  problem  since  a  region  will  not  easily 
grow  beyond  an  existing  jump  edge.  On  the  other  hand,  crease  edges  are  not  easily 
accounted  for  since  surface  fitting  alone  may  grow  beyond  simple  discontinuities  in 
orientation.  Therefore  code  has  b»  ?n  added  that  can  be  used  to  limit  region  growth 
based  upon  a  mask.  Such  a  mask  would  be  made  based  upon  crease  and  jump  edges 
in  the  image.  This  option  has  been  tested  and  works  functionally.  It  could  prove 
very  effective  when  integrating  other  types  of  sensor  information.  For  example,  edge 
information  provided  by  a  passive  IR  image  could  be  taken  into  account. 

So  far  we  have  seen  that  the  target  regions  can  be  segmented  into  separate  planar 
surface  regions,  and  the  reconstruction  image  visually  resembles  the  original  image  with 
much  less  noise.  But  we  have  found  that  oveisegmentation  often  occurs,  with  too  many 
subregions  being  formed.  Post  processing  to  merge  similar  regions  that  are  adjacent  would 
help  eliminate  this  problem,  Crease  edge  detection  will  be  used  at  this  point  to  help 
determine  if  two  regions  are  distinct  or  if  they  need  to  be  merged, 

2.3  Robust  Planar  Surface  Estimation 

This  investigation  was  initiated  when  it  was  observed  that  the  planar  fits  that  were  being 
found  for  sections  of  range  imagery  were  highly  erratic.  Error  in  surface  estimates  may  be 
due  to  several  different  causes. 

e  Expected  variance. 

•  Incorrect  modelling  of  the  noise  process. 

•  Data  samples  that  do  not  come  from  a  single  underlying  distribution. 


«  Invalid  parametric  model, 

It  is  always  desirable  to  chose  a  technique  which  has  low  expected  variance  but  other 
points  must  be  considered  as  well.  An  estimation  technique  which  v^rks  well  under  specific 
conditions  but  degrades  rapidly  with  very  minor  changes  in  the  ura  arlying  distribution  is 
not  useful  when  there  is  some  doubt  as  to  the  form  of  the  noise  process.  Parametric  surface 
estimation  techniques  generally  assume  all  the  data  points  belong  to  the  same  underlying 
distribution.  For  instance  if  g&ussian  noise  is  the  only  problem,  then  the  generally  accepted 
estimation  technique  is  least  squares  regression.  When  corrupting  data  which  belongs  to  a 
different  underlying  distribution  is  present  then  the  original  surface  estimation  technique 
will  have  degraded  performance.  A  problem  also  arises  when  it  is  suspected  that  the 
parametric  model  is  inaccurate.  For  instance,  does  all  the  data  really  come  from  one 
plane,  or  does  part  of  it  belong  to  another  surface.  Some  estimation  techniques  can  handle 
the  addition  of  corrupting  data  more  readily  than  others.  Since  some  LADAR  returns 
may  be  corrupted  and  slight  missegmentations  may  also  occur,  it  is  important  to  find  a 
technique  that  is  robust  under  these  conditions. 

In  this  section  we  discuss  briefly  the  most  promising  techniques  for  planar  estimation  in 
the  presence  of  considerable  noise  and  some  corrupting  data.  Comparisons  of  their  behavior 
on  controlled  data  are  also  detailed.  These  methods  include:  least  squares  regression,  M- 
estimate  regression,  and  the  removal  of  outliers  based  upon  a  weighting  function.  Each  of 
these  techniques  is  described  and  compared  in  the  following  sections. 

Techniques 

Initially  the  data  is  assumed  to  come  from  a  plane  with  additive  noise.  The  position 
variable  is  x,  --  (I,x,,ta),  and  the  corresponding  range  observation  is  z,.  Then  we  say 

Zi  -  •  /?“  i  x,/?"  +  y,7?j  i  e,  —  x,/?°  I  c, 

where  the  e,  are  independent,  identically  distributed  noise  processes  with  mean  0  and 
standard  deviation  cr,  and  ft0  (ft®,  fit*  P2)  are  the  planar  coefficents.  If  there  are  n  data 
points  then  we  have 

z  X(  f  f  e 

where  X  is  the  matrix  of  x, 


Least  Squares:  The  goal  in  least  squares  is  to  estimate  ft  so  as  to  minimize  the  difference 
between  the  predic  ted  range  values,  x,ft,  and  the  observed  values,  z,,  i.e. 

ri 

min  N  ( z,  x,/f) “ 

ft  ( X  1  X )  1  X  r  r. 

f> 


This  minimum  occurs  when 


As  one  would  hope  the  expected  value  of  /?  is  /3°.  The  covariance  of  is 

Cov(0)  =  o*(XTX)~l. 

Thus,  if  the  region  to  be  fit  is  a  square  of  area  n  centered  at  the  origin,  then 

Var(^o)  =  o2(n 
Var(ft)  « 

n 

Var(^j)  »  <?*(—) 
n. 

So  if  a  is  large,  as  is  often  the  case  in  range  imagery,  then  a  large  number  of  pixels  is  needed 
in  each  region  to  obtain  an  accurate  fit.  The  variance  of  (3Q  will  depend  upon  the  location 
of  the  origin  relative  to  the  center  of  the  region  being  fit,  with  a  larger  distance  inducing  a 
larger  variance.  The  slope  variances  do  not  change  with  translation,  but  they  will  change 
as  the  shape  of  the  region  is  altered.  For  example,  a  region  which  is  long  in  the  x  direction 
but  narrow  in  the  y  direction  will  have  a  small  variance  for  pi  but  a  larger  variance  for  /?2. 
A  more  compact  region  that  has  the  same  area  would  have  a  larger  variance  for  and  a 
smaller  variance  for  j3j. 

However,  sometimes  the  observed  data  does  not  exhibit  the  aforementioned  behavior. 
In  other  words  the  above  description  of  e,-  is  incorrect.  Specifically,  certain  data  points 
may  be  corrupted  by  more  than  the  expected  noise.  Since  least  squares  tries  to  minimize 
the  squared  error,  undue  weight  may  be  given  to  points  that  are  out  of  line  with  the  rest 
of  the  data.  For  example,  consider  the  least  squares  regressions  displayed  in  Figure  I. 


Figure  1:  Least  squares  linear  regression  when  some  of  the  data  is  corrupted. 


There  are  a  variety  of  ways  of  dealing  with  this  type  of  problem.  One  method  which 
retains  all  the  data  is  based  upon  M-estiinates. 

Regression  M-estimate:  In  least  squares  regression  the  objective  is  to  minimize  the 
sum  of  the  squared  errors.  This  gives  great  importance  and  control  to  those  data  points 
that  normally  would  have  a  large  error.  If  the  sum  of  the  squares  is  replaced  by  a  less 
rapidly  increasing  function  of  the  residuals  an  M-estimate  is  produced.  The  problem  then 
becomes 

min  )  />(z,  x,/i) 

i)  - 
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for  some  coat  function  p(-).  The  minimum  can  be  characterized  by  the  following  equations: 


where 


^x{^£) 


0,  k-  0,1,2 
0, 


V>(*)  =  />'(*), 

X  (x)  =  xt/>(x)  -  p(x) 


For  the  test  case  used  here,  the  cost  function  rwitches  from  square  to  linear  once  the 
residual  becomes  too  large.  In  particular 

^’(x)  =  max(--e,  min(c,x))  with  c  =  1.5. 

The  constant  c  was  chosen  based  upon  the  experience  of  experts  in  linear  surface  esti¬ 
mation  Additional  testing  of  this  technique  on  LADAR  data  is  needed  for  an  optimal 
setting.  Since  this  is  a  nonlinear  function,  the  regression  must  be  solved  by  iteration.  The 
algorithm  used  is  detailed  in  Huber  |5j.  For  more  specifics  on  the  asymptotic,  behavior  of 
M-estimates  in  regression  see  Huber  [6]  as  well, 

Another  technique  for  dealing  with  data  points  that  are  suspected  of  being  corrupted 
is  to  remove  them  entirely  from  the  least  squares  computation.  First  though,  these  data 
points  must  be  identified  as  outliers. 


Outlier  identification:  The  identification  of  outliers  in  a  sample  usually  proceeds  via  a 
variety  of  tests  with  a  good  deal  of  human  subjectivity  involved.  Most  methods  start  with 
a  least  squares  regression  using  all  the  points.  Then  the  residuals  from  that  regression  and 
knowledge  of  the  design  matrix  X  are  used  to  find  outliers. 

The  residuals  from  least  squares, 

e  =  -  Xp  (/  X(XTX)~lXT)z  -----  G’z, 

have  a  variance  that  changes  according  to  the  location  of  x,.  For  example,  the  addition  of 
a  constant  ffsot  to  all  '.he  data  samples  will  affect  the  residuals  e.  For  this  reason  it  is  not 
wise  to  compare  the  residuals  directly  with  each  other.  Instead  there  is  a  “normalised* 
version  of  the  residue!, 

c. 

r8  =  — 7= 

OyJg,t 

which  has  a  "arianoe  of  1  so  that  the  values  are  comparable.  If  the  variance  of  tiie  noise, 
a2,  is  not  known,  then  it  can  be  estimated  from  the  data  a- 


n  3 
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A  large  normalized  residual  is  not  the  only  indicator  of  a  bad  data  point.  In  the  first, 
plot  in  Figure  1  the  residual  at  the  last  point  is  not  very  large  but  it  is  still  an  outlier.  This 
point  has  undue  influence  on  the  results  of  the  regression  and  is  called  a  leverage  point. 
The  values  on  the  diagonal  of  the  matrix  G,  ga ,  are  related  to  the  leverage  at  point  *'.  As 
ga  approaches  0,  the  leverage  increases. 

There  are  several  ways  of  combining  the  values  of  r,  and  gie-  into  an  influence  function 
that  measures  the  effect  of  point  i  on  the  total  regression  as  compared  to  the  other  points. 
Two  of  them  are: 

(1)  Normalized  residual  #i(rd,  £h«)  =~  rf ,  and 

(2)  Cook’s  function  jf,«)  =  -(-1 — — ) 

3  ga 

Other  examples  of  influence  functions  are  given  in  [3]. 

Those  data  points  which  have  abnormally  large  values  of  the  influence  function  are 
removed  from  the  data  set  and  least  squares  regression  is  repeated  without  them.  This 
process  can  be  repeated  as  needed.  Several  questions  arise  when  using  this  strategy: 

•  What  is  meant  by  “abnormally  large”? 

•  How  many  points  can  be  removed  at  once? 

•  When  should  the  process  be  terminated? 

The  variance  of  the  normalized  residual  is  one  and  the  mean  is  zero,  so  if  the  errors  are 
all  gaussian  we  will  have  |r,|  <  2.57  with  probability  .99.  Therefore  a  value  of  Hi  greater 
than  6.6  is  highly  suspect.  Cook’s  function  does  not  have  a  p re-determined  variance  which 
is  constant  over  all  the  data.  It  measures  the  effective  change  in  /?  that  would  occur  if 
the  s  -th  data  point  were  removed.  The  removal  of  most  points  will  have  little  effect  on 
the  regression  estimate.  Those  points  whose  removal  causes  a  change  that  is  significantly 
larger  than  average  are  potential  outliers.  These  points  are  identified  by  estimating  the 
variance,  s2,  of  r,y(  1  —  gu)/ga  over  the  sample  data.  When  H2  >  6.6s2  there  is  probably 
an  outlier  present. 

for  the  highest  accuracy,  it  would  be  best  to  remove  only  one  point  -luring  each  pass. 
This  allows  the  residuals  to  readjust  and  avoids  forcing  the  estimate  to  remain  the  same 
because  all  non-supporting  data  has  been  removed  immediately.  This  is  a  veiy  time  con¬ 
suming  process,  however,  so  there  is  a  tradeoff  between  time  and  the  potential  for  ruining 
the  computation.  The  potential  for  eliminating  a  large  set  of  correct  points  is  avoided  if 
the  number  of  points  that  can  be  removed  per  pass  is  fixed  to  a  certain  percentage  of  the 
points  currently  present.  In  these  tests  at  most  10%  of  the  points  are  targeted  for  removal 
in  the  each  pass. 

The  outlier  removal  can  continue  until  no  more  points  meet  the  outlier  criterion.  This 
is  a  dangerous  practice,  however,  since  the  number  of  points  can  be  become  so  small  that 
the  results  are  meaningless.  It  is  important  therefore  to  limit  the  total  number  of  points 
that  can  be  eliminated,  as  well  as  the  number  of  outliers  removed  at  each  stage.  For  these 
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tests,  the  limit  on  total  number  of  points  removed  is  based  upon  the  size  of  the  perimeter 
—  since  it  will  often  be  the  perimeter  points  which  are  suspect 

Tests 

In  order  to  teat  the  planar  estimation  routines  described  above,  several  sets  of  synthetic 
data  were  generated.  Parameters  were  varied  to  produce  the  data  sets.  These  parameters 
included: 

•  Image  Size:  Either  5x5,  10x10,  or  15x15. 

•  Expected  Noise:  a  —  1,  or  5. 

•  Extra  Corruption:  0%  or  10%  of  the  points  had  added  noise  with  standard  deviation 
a  =  20.  The  corruption  was  either  Gaussian  or  a  one-sided  Gaussian. 

•  Angle  of  Plane:  Either  perpendicular  to  the  line  of  sight  0°  or  at  45°. 

The  specifics  of  the  different  data  sets  are  contained  in  Table  1. 

Each  of  the  fitting  techniques  described  in  the  previous  sections  was  run  on  all  the 
synthetic  data.  Table  2  shows  the  labels  that  are  used  to  denote  the  different  algorithms. 
To  score  each  run,  the  interior  angle  between  the  computed  normal  and  the  true  normal 
was  calculated.  The  closer  this  angle,  6,  is  to  zero  the  better  the  result.  The  scores  from 
each  of  the  runs  are  shown  in  Tables  3  and  4.  The  different  algorithms  were  then  compared 
by  using  the  Wilcoxon  signed-ranks  pairwise  hypothesis  test  (see  PM  II  monthly  report 
August  1986). 

It  should  be  noted  that  other  tests  that  were  considered,  but  were  not  performed  due  to 
lack  of  time.  For  instance,  a  test  of  data  which  had  corrupted  pixels  along  one  boundary, 
simulating  the  effect  of  slight  missegmentations,  could  have  been  performed. 

Conclusions 

First  it  should  be  noted  that  when  the  noise  was  increased  to  a  —  5,  as  in  data  set  DS4, 
the  planar  fits  r  a  5x5  region  were  so  distorted  that  they  were  useless.  Because  of  this, 
future  tests  with  a  =  5  were  not  performed  on  the  5x5  images  and  the  results  from  the 
5x5  images  in  DS4  were  ignored.  Also  a  check  of  the  estimated  means  and  variances  of 
9  between  data  sets  DSl,  DS2,  and  DS3,  and  their  counterparts  in  DS7,  DS8,  and  DS9 
shows  that  the  viewing  angle  of  the  plane  does  not  significantly  alter  the  behavior  of  the 
algorithms.  The  slight  improvement  in  the  angular  results  when  the  viewing  angle  is  45° 
is  due  to  slope  to  angle  transformation  rather  than  any  actual  performance  improvement. 
The  number  of  pixels  which  can  be  seen  on  the  plane  is  a  much  more  important  factor 
in  determining  algorithm  performance.  Generally  of  course,  for  a  fixed  size  surface,  the 
number  of  viewed  pixels  will  be  proportional  to  the  viewing  angle. 

Table  5  presents  the  results  from  the  pairwise  hypothesis  tests,  and  the  conclusions 
are  summarized  here.  Each  pairwise  test  that  is  performed  has  a  potential  “winner”  and 


*lo«er" .  When  th*  significance  level  becomes  small  enough,  there  is  evidence  to  support  this 
winner- loser  division.  If  the  level  is  too  large,  neither  algorithm  can  be  called  a  “winner”. 
In  fact  a  result  as  high  as  0.3  is  meaningless,  and  comparisons  that  produced  such  results 
are  not  presented  in  the  table.  On  the  other  hand,  the  appearance  of  an  algorithm  pair  in 
the  table  does  not  guarantee  that  one  is  better  than  the  other.  Confidence  in  this  type  of 
conclusion  must  be  based  upon  the  significance  level.  For  instance,  in  data  set  DSl  we  sec 
that  the  M -estimate  (M)  performed  better  than  Normalised  residual  outlier  removal  (N) 
and  Least  Squares  (LS)  with  a  fair  amount  of  certainty,  and  better  than  Cook’s  function 
outlier  removal  (C)  with  less  certainty.  But,  no  one  of  N,  LS,  or  C  can  be  said  to  have 
performed  better  than  any  of  the  others. 

A  brief  scan  of  the  results  of  the  Wilcoxon  hypothesis  tests  on  the  separate  data  sets, 
shows  that  the  M-estimate  is  often  the  top  performer,  while  least  squares  is  often  the  worst. 
As  the  noise  increases  the  least  squares  estimate  appears  to  perfoxm  relatively  better  than 
at  low  noise  levels.  In  fact  when  only  the  expected  noise  is  present  least  squares  is  the 
best  performer,  but  it  is  still  at  a  disadvantage  when  there  is  extra  corruption  in  the 
image  beyond  the  expected  noise.  Since  least  squares  is  the  maximum  likelihood  estimate 
when  the  noise  distribution  is  truly  Gaussian  it  must  be  concluded  that  quantization 
has  significantly  distorted  the  distribution.  This  is  hardly  surprising  in  the  case  where 
o=l.  The  performance  of  the  M-estimate  depends  upon  the  value  of  the  threshold,  c, 
which  remained  fixed  in  these  teats.  Ideally  this  threshold  would  vary  based  upon  the 
characteristics  of  the  input  data.  (There  axe  allusions  to  this  in  the  literature,  but  we  have 
not  yet  seen  any  specifics.)  Therefore,  the  M-estimate  may  perform  even  better  than  is 
demonstrated  in  these  experiments. 

When  all  the  data  sets  arc  merged  together  a  very  definite  ordering  of  algorithm  per¬ 
formance  appears. 

1.  M-estimate. 

2.  Normalized  residual  outlier  removal. 

3.  Cook’s  function  outlier  removal. 

4.  Least,  squares. 

When  using  parametric  surface  estimation  alone,  M- estimates  is  the  most  robust  of  the 
four  techniques  tested.  M-estimates  is  iterative  and  can  be  time  consuming,  so  it  might 
be  possible  to  get  similar  results  by  (i)  using  Gaussian  smoothing  to  remove  part  of  the 
quantization  error,  (ii)  applying  a  replacement  type  filter  to  remove  extreme  points,  and 
(iii)  proceeding  with  least  squares. 


Table  1:  Data  sets  used  for  planar  estimation  tests 
(See  text  for  options.) 


Data 

Angle 

Expected 

Gaussian 

No. 

of  Images/Size  | 

Set 

NoIsj 

Corruption 

5x5 

10x10 

15X15 

DSl 

0° 

a  =  1 

5 

S 

S 

DS2 

a  --  1 

10% 

5 

5 

S 

DS3 

0° 

<7=1 

10%  one-sided 

5 

S 

5 

DS4 

0° 

BSE 

5 

5 

5 

DSS 

0° 

cs5 

10% 

5 

5 

OS6 

0° 

ESI 

10%  one-sided 

n 

S 

DS7 

45° 

<7=1 

5 

5 

OS  7 

45° 

a  =  1 

10% 

5 

5 

5 

OS7 

45° 

<7  =  1 

10%  one-sided 

5 

9 

5 

Table  2:  Labels  of  the  plane  fitting  algorithms. 


LS:  Standard  least  squares. 

M:  Robust  M-esilmate  regression. 

N:  Outlier  removal  bared  on  normalized  residuals. 

C:  Outlier  removal  bared  on  Cook’s  function. 
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Table  3:  Angie,  0  (in  degrees)  between  result  normal  and  actual  normal 
The  mean,  m,  and  estimated  variance,  sd,  of  the  angle  are  given  for  each  size  grouping. 


IS 

M 

N 

MB 

■m 

LS 

M 

N 

C 

LS 

M 

N 

c  I 

5.71 

4.50 

5.55 

5.71 

63.43 

24.28 

19.21 

21.58 

Ksa 

ll.ftS 

13.12 

13.24 

13.24 

56.24 

26.87 

18.27 

17.67 

BS 

4.57 

1.57 

4.57 

4.57 

34.59 

0.19 

22.90 

19.51 

22.99 

23.40 

25.24 

25.24 

18.53 

15.8S 

10.57 

10.57 

41.51 

13.34 

0.50 

9. SC 

82.40 

19.03 

16.98 

19.98 

12.7S 

0.10 

12.75 

12.75 

31.31 

14.14 

14.76 

14.76 

49.75 

16.31 

24.01 

31.13 

10.02 

4.7 

4.85 

5.04 

10.94 

5.54 

10.07 

5.80 

45.48 

13.86 

15.78 

10.57 

14.93 

S.07 

16.01 

4,70 

42.40 

16.18 

15.39 

8.80 

19.69 

9.67 

19.88 

8.62 

4.80 

S.31 

3.81 

2.06 

4.03 

S.25 

4.70 

4.03 

0.85 

1.86 

1.88 

2.24 

2.23 

1.18 

1.45 

1.80 

1.23 

1.66 

11.61 

4.69 

4.57 

5.52 

1.94 

0.04 

3.03 

2.49 

0S2 

5.93 

1.24 

2.47 

OSS 

0.00 

5.62 

6.08 

4.2S 

2.87 

2.95 

2.67 

2.57 

37.45 

0.95 

21.50 

35.23 

22.42 

1.15 

1.63 

1.63 

2. S3 

2.72 

1.04 

1.78 

19.46 

1.37 

3.30 

3.78 

12.26 

3.59 

1.40 

1.65 

2.78 

1.0« 

2.25 

1.30 

2.35 

1.20 

2.07 

0.47 

24.5 

8.78 

2.80 

2.07 

6  10 
8.67 

0.57 

14.39 

11.66 

6.85 

3.18 

2.12 

3.13 

2.08 

2.08 

1  80 

0.3S 

0.60 

0.32 

1.20 

1.20 

1.14 

5.96 

0.79 

0.98 

1.01 

1.70 

0.7* 

1.31 

0.52 

0.71 

1.17 

11.03 

1.80 

7.53 

7.99 

1.40 

1.04 
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4.59 

0.54 

1.07 

0.85 

7.71 

2.40 

1.61 

3.32 

o.«e 

0.30 

0.10 

6.44 

2.80 

4.61 

5.77 

4.64 

3.62 

2.12 

3.73 

0.30 

0.27 

0.72 

8.00 

2.70 

1.05 

1.57 

7.04 

0.34 

2.65 

4.82 

0.80 

0.63 

0.6b 

0.33 

0.50 

0.28 

0.71 

0.62 

4.31 

3.56 

1.56 

1.12 

1.73 

1.62 

2.10 

2.07 

7.46 

2.41 

1.79 

1.31 

2.98 

2.62 

4.17 

2.55 

LS 

Ml 

N 

C 

LS 

M 

N 

C 

LS 

M 

N 

C 

63.07 

62.66 

62.27 

63.07 

34.58 

44.41 

34.58 

34.58 

52.50 

56.53 

52.58 

45.53 

32.53 

36  83 

43.76 

43.83 

53.17 

48.42 

53.17 

59.60 

47.18 

13.14 

48.77 

10.12 

48.28 

10.50 

49.36 

11.30 

15.55 

20.77 

21.41 

29.03 

15.10 

8.74 

3  86 

4.70 

11.01 

15.19 

27.66 

24.80 

6  54 

13.23 

3.52 

12.67 

13.22 

22.70 

26.79 

23.37 

15.92 

2.08 

9.41 

10.58 

21.02 

23.74 

34.89 

33.61 

DS5 

11.43 

4.50 

13.86 

10.93 

DS6 

21.70 

23.20 

36.24 

33.06 

10.07 

14.67 

12.54 

12.01 

15.84 

8.08 

6.61 

2.49 

25.21 

18.82 

15.25 

7.93 

15.80 

23.68 

34.25 

24.94 

11.86 

8.16 

18.35 

26.41 

9.50 

11.63 

13.35 

13.30 

13.80 

5.61 

19.22 

4.98 

21.32 

13.65 

22.45 

9.73 

14.49 

2.85 

10.64 

6.09 

14.11 

9.35 

13.58 

10.83 

16.67 

6.75 

14.30 
7. 78 

20.38 

11.19 

17  93 
10  63 

5.98 

4  32 

5.50 

8.34 

13  27 

7.05 

8.67 

12.83 

9.47 

6.06 

7.19 

7  OS 

2.72 

4.87 

4.08 

5.15 

7.33 

8.15 

2.73 

4  19 

3.52 

5.26 

6.77 

7  61 

552 

5  24 

6  54 

5  31 

4.23 

3.55 

13  58 

9.33 

6  66 

3.76 

7.14 

10  4  S 

m 

id 


10 

x 

10 


m 

Id 

15 

x 

15 


DS1 


m 

sd 


m 

sd 


10 

x 

10 


DS4 


m 

sd 


15 
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Table  4:  Angle,  8,  (in  degrees)  between  result  normal  and  actual  normal  (continued). 
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Table  5:  Results  of  the  Wilcoxon  signed-ranks  hypothesis  tests.  Each  triple  consists  of: 
the  “winner”,  the  “loser”,  and  the  significance  level.  Smaller  significance  levels  indicate 
higher  confidence  in  the  “win- lose”  decision.  Triples  with  levels  >  0.3  are  not  shown. 
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2.4  Ground  Plane  Removal 

This  thocholda  used  for  the  region  growing  segmentation  are  computed  directly  from  the 
image,  but  experience  has  shown  that  they  are  set  inappropriately  when  a  large  amount  of 
background  is  Included  in  the  computation.  This  is  to  be  expected  since  the  background 
is  quite  often  not  as  smooth  and  regular  as  the  target  itself.  Therefore,  there  is  a  need 
to  eliminate  most  of  the  background  in  an  image  before  computing  the  surface  fitting 
thresholds. 

The  goal  of  this  effort  was  not  to  totally  segment  the  target  object  from  the  background, 
but  to  identify  areas  that  can  be  reasonably  expected  to  belong  to  the  background.  Two 
major  techniques  were  investigated  in  this  regard. 

•  Create  closed  boundaries  from  edge  pixels. 

•  Model  the  background  as  a  plane  and  strive  to  find  this  plane. 

The  first  technique  used  edge  information  along  with  edge  linking  to  try  to  produce 
closed  boundaries.  This  was  found  to  be  unreliable  since  linking  v/aa  rarely  complete 
and  the  lower  edge  between  target  and  ground  plane  was  often  not  detected.  The  linking 
techniques  based  upon  edge  magnitude  and  direction  may  prove  useful  though  when  trying 
to  identify  crease  edges.  Figure  2  shows  an  example  of  the  edge  detection  on  synthetic 
noise  free  data.  The  failure  to  detect  the  bottom  edge  of  the  target  occurs  even  in  this 
noise  free  environment. 

The  second  technique,  based  on  a  background  model  rather  than  edges,  was  more 
effective.  One  of  the  largest  major  components  cf  the  background  is  the  ground  plane. 
Therefore  efforts  were  directed  at  finding  and  identifying  a  planar  region  that  occupies  a 
major  portion  of  the  boundary  of  the  image,  and  extrapolating  this  information  to  the 
interior  of  the  image. 

Two  different  approaches  to  finding  and  removing  the  ground  plane  were  investigated. 
The  first  involved  simple  accumulation  of  evidence.  It  was  assumed  that  the  majority  of 
the  boundary  of  the  image  is  part  of  the  background.  Therefore  a  planar  fit  was  computed 
for  a  boundary  point,  and  then  any  point  within  the  image  that  lay  close  to  this  plane 
accumulated  evidence  that  it  was  also  a  part  of  the  ground  plane.  Once  this  was  done  for 
all  the  boundary  points  in  turn,  those  points  with  the  highest  accumulation  of  evidence 
would  be  background  points.  Unfortunately,  this  technique  was  so  time  consuming  relative 
to  the  surface  extraction  itself  that  it  was  judged  infeasible. 


Figure  2:  Edge  detection  on  a  synthetic  range  image  of  a  T-55  tank. 
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The  second  Approach  to  removing  the  ground  plane  is  based  on  similar  assumptions, 
but  usee  clustering  techniques  instead  of  accumulation.  Briefly,  the  planar  normal  and 
intercept  are  approximated  at  each  point  in  the  image  and  the  normals  are  clustered  to 
find  parallel  planar  regions.  The  cluster  moat  likely  to  contain  the  background  plane  is  then 
separated  based  upon  intercept  value  so  that  object  regions  parallel  to  the  backgound  plane 
are  not  removed.  This  approach  is  much  faster  than  the  previous  one,  and  is  described  in 
more  detail  in  what  follows. 

The  normal  at  each  point  is  represented  by  the  x.  and  y  partial  derivatives  and  the 
intercept  is  computed  (mod  256)  using  the  center  of  the  image  as  the  origin.  The  normals 
are  then  clustered  in  two-dimensional  space  using  a  convergent  k-meana  technique.  Since 
we  want  to  separate  out  one  particular  area  in  2-D  space,  a  5-means  clustering  was  used. 
This  allows  for  clusters  on  each  side  of  the  principal  cluster.  Brief  tests  with  other  numbers 
of  means  have  all  produced  inferior  results.  The  clustering  routine  currently  available 
employs  the  Euclidean  norm,  but  we  would  expect  improved  results  if  the  Mahalanobis 
norm  were  substituted. 

In  order  to  determine  which  of  the  clusters  contains  the  ground  plane  the  assumption 
of  backgound  occupying  the  majority  of  the  boundary  is  again  used.  Other  types  of 
information  such  as  the  depression  angle  or  terrain  maps  could  be  used  to  enhance  this 
determination.  Once  this  clurter  is  identified  it  is  further  divided  based  upon  the  intercept 
values.  The  intercepts  are  histogrammed  to  find  the  major  background  intercept  region, 
so  that  the  other  parallel  areas  can  be  removed  from  consideration.  This  is  done  by 
finding  the  maximum  histogram  value  and  moving  outward  until  the  value  drops  below 
a  given  percentatge  of  the  maximum.  Points  with  the  appropriate  normal  and  intercept 
value  which  are  connected  to  the  boundary  are  identified  as  potential  background  points. 
Finally  those  areas  not  currently  identified  as  backgound  are  classified  as  either  an  area  of 
interest  based  upon  some  external  criterion  or  added  to  the  background  region. 

Figures  3  and  4  show  various  stages  of  the  ground  plane  removal  performed  on  synthetic 
imagery.  The  original  image  in  Figure  3  is  a  T-55  tank  range  image  generated  from  ERIM 
wire  frame  model  with  ground  plane  added  at  a  depression  angle  of  10°.  Gaussian  noise 
was  added  to  the  image  to  create  SNR  =  4.  The  sequence  shows  the  5  different  normal 
clusters,  the  chosen  background  cluster  with  points  removed  based  upon  intercept  value, 
the  connected  portion  of  the  background,  and  finally  the  foreground  area  which  will  be 
handed  off  for  segmentation. 

The  original  image  in  Figure  4  is  an  M-109  viewed  at  a  depression  angle  of  10°  and  the 
corresponding  ground  plane.  No  noise  was  added  to  this  image.  In  this  case,  the  sequence 
shows  the  background  cluster  both  before  and  after  separation  based  upon  the  intercept, 
and  the  mask  to  be  used  for  segmentation.  The  only  substantial  differences  between  the 
final  mask  and  the  original  object  are  the  artifacts  due  to  the  range  ambiguities.  The  mask 
is  slightly  larger  than  the  real  object  because  of  averaging  effects. 

To  give  an  idea  of  the  power  of  the  surface  segmentation,  Figure  5  shows  the  segmen¬ 
tation  of  the  original  noise  free  T-55  image  and  the  M-109  image.  The  sawtooth  edge 
occuring  in  the  M-109  segmentation  is  due  in  part  to  quantization  effects.  The  gun  barrel 
of  the  T-55  is  included  as  part  of  one  of  the  turret  planes  since  the  narrow  barrel  produces 
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(8)  (H)  (i) 


Figure  3:  Ground  Plane  Removal:  (a)  original  image:  T-55  at  SNR  =  4,  gun  barrel 
clipped  to  reduce  image  size,  (b-f)  points  included  in  each  of  the  5  normal  clusters,  (g) 
points  in  ground  plane  duster  minus  those  with  the  wrong  intercept  values,  (h)  connected 
component  of  ground  plane,  (i)  final  mask  identifying  foreground. 
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Figure  4:  Ground  Plane  Removal:  (a)  original  image:  M-109  without  noise,  (b)  ground 
plane  cluster,  (c)  after  removal  of  points  based  upon  intercept,  (d)  final  mask  of  foreground. 


Figure  5:  Surface  Segmentation  Examples  without  noise:  (a)  T-55,  (b)  M-109 


degenerate  planar  fits  by  itself.  Thb  might  be  corrected  by  also  allowing  quadratic  fits. 


2*5  Smoothing  Prior  to  Computing  Derivatives 

First  and  second  order  partial  derivatives  are  often  used  in  range  image  analysis,  but  mea¬ 
suring  them  reliably  is  difficult.  Edge  detection  may  be  viewed  as  the  characterization  of 
intensity  changes  within  an  image.  Thus,  effective  edge  detection  necessitates  the  reliable 
estimation  of  the  differential  behavior  of  the  image  function.  Specifically,  derivatives  of 
various  types,  order  and  scale  have  to  be  accurately  calculated  for  edge  detection  to  func¬ 
tion  correctly.  Surface  normals  and  Gaussian  and  mean  curvature  are  also  functions  of 
partial  derivatives  and  thus  have  the  same  reliability  problems. 

Numerical  differentiation  of  a  sensed  image  function  is  an  ill-poeed  problem  because  its 
solution  does  not  depend  continuously  on  the  data.  This  may  be  directly  attributed  to  the 
random  perturbations  introduced  by  the  previous  sensing  and  sampling  processes.  Thus, 
before  edge  detection  can  be  reliably  performed,  the  differentiation  problem  must  become 
well-posed  and  stable. 

As  argued  in  [11]  ill-posed  problems  can  be  be  solved  using  the  regularization  method¬ 
ology  originally  suggested  in  [9,15].  Briefly,  the  general  approach  is  to  restrict  the  space  of 
possible  solutions  to  subspace  by  imposing  supplementary  constraints.  Of  the  many  possi¬ 
ble  regularization  techniques,  one  that  is  natural  for  edge  detection  on  range  images  is  that 
of  Reinsch  [12]  considered  for  numerical  differentiation.  In  the  context  of  edge  detection, 
this  approach  suggests  optimal  filters  and  their  approximations  that  can  be  applied  to  the 
noisy  data  prior  to  computing  derivatives  (ll). 

In  the  literature,  various  criteria  have  been  proposed  for  directing  the  determination  of 
“optimal”  regularization  filters.  Basically,  these  criterion  center  upon  a  tradeoff  between 
the  degree  of  regularization  of  the  filtered  data  versus  the  closeness  of  filtered  data  to 
the  original  data.  Once  the  criterion  has  been  selected,  a  regularization  filter  must  be 
determined  for  each  specific  type  of  derivative. 

For  jump  edges,  a  Gaussian  operator  is  a  good  approximation  [ll]  to  the  optimal  filter 
derived  using  a  calculus  of  variations  technique.  Although  the  regularization  approach 
suggests  the  near  optimality  of  Gaussian  filters,  the  variance  of  the  Gaussian  function 
corresponding  to  the  scale  of  the  filter  must  be  specified.  At  least  two  or  three  methods 
for  determining  the  scale  have  been  suggested  in  the  computer  vision  literature  [11,18)  but 
we  are  not  aware  of  any  that  have  been  applied  to  infrared  or  range  images. 

In  addition  to  finding  the  optimal  scale  for  a  Gaussian  filter  there  is  the  option  of 
applying  the  filter  at  various  scales.  In  an  image,  changes  of  intensity  take  place  at  many 
spatial  scales  depending  on  their  physical  origin.  Consequently,  a  multiscale  analysis, 
tracing  the  behavior  of  some  feat  ure  of  the  signal  across  scales,  can  reveal  useful  information 
about  the  nature  of  the  underlying  physical  process.  For  example  [19],  spatial  coincidence 
at  all  scales  of  zero  crossings  in  the  Laplacian  of  Gaussian  filtered  images  may  confirm 
a  physical  “edge”  distinct  from  surface  markings.  It  should  be  emphasized  that  ic  is  not 
only  necessary  to  detect  and  describe  changes  in  a  range  image  at  different  scales,  but  in 
addition,  much  useful  information  can  be  obtained  by  combining  descriptions  across  scales. 
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Currently  a  Gaussian  filter  is  being  applied  to  the  range  data  prior  to  computing  the- 
derivatives  needed  for  Gaussian  and  mean  curvature.  The  results  of  this  filtering  are 
particularly  obvious  on  synthetic  data.  Quantization  effects  can  easily  be  seen  on  the 
HK-Sign  map  for  a  plane  which  is  at  an  angle  to  the  line  of  sight.  Instead  of  labelling 
the  entire  plane  as  a  flat,  periodic  stripes  occur  across  the  plane  of  other  surface  types. 
Regularization  filtering  helps  to  minimize  this  quantization  effect  as  well  as  suppressing 
the  noise. 

2.6  Crease  Edge  Detection 

During  surface  segmentation  it  is  possible  that  a  single  surface  becomes  subdivided  into 
several  separate  regions.  If  these  separate  regions  are  used  when  forming  the  object  graph 
then  any  matching  will  be  unlikely  to  succeed.  Therefore  a  method  must  be  found  for 
deciding  when  two  regions  belong  to  the  same  surface  as  opposed  to  two  separate  surfaces. 

We  are  proposing  to  make  decisions  on  the  distinctness  of  two  regions  based  upon  the 
presence  or  absence  of  crease  edges.  If  a  crease  edge  is  present  then  the  two  regions  remain 
separate.  Otherwise  the  two  adjacent  regions  are  merged  and  considered  as  one  surface. 

Jump  edges  tend  to  be  easily  identified  by  the  surface  growing  algorithm,  but  crease 
edges  are  harder  to  find.  Mitiche  and  Aggarwal  [8]  have  pt  opened  a  method  to  crease 
edge  detection  which  has  the  advantage  of  being  based  upon  a  known  noise  model.  This 
technique  has  been  implemented  and  is  in  the  process  of  being  tested  and  integrated  with 
the  surface  segmentation  software. 

This  approach  is  based  on  fitting  planes  in  various  directions  in  a  local  neighborhood 
and  computing  the  angles  between  the  planes.  Let  g{xQ,y0)  be  the  depth  map  value  at 
(®o,J/o)  and  N  be  an  appropriate  sized  neighborhood  around  (xo,  J/o)-  Let  N  be  divided 
into  two  contiguous  planar  surfaces,  ct  and  c2,  n j,  and  n2  being  the  surface  normals  of 
the  planes  best  fit  (say  in  least  square  sense)  to  the  points  in  c2  and  c2  respectively.  The 
collection  of  points  in  N  together  with  these  best  fit  plants  will  be  referred  to  as  a  partition 
of  N  and  is  denoted  as  R.  Suppose  n  is  the  number  of  directions  in  which  crease  edges 
have  to  be  determined.  Iri  our  implementation  n  is  4,  corresponding  to  horizontal,  vertical, 
left  diagonal,  and  right  diagonal  directions.  For  each  direction,  planes  are  fit  according 
to  some  selected  criterion  to  the  corresponding  two  regions  in  N  thus  determining  the 
partitions  R+,  i  —  1, . . .  ,  n.  The  angles  between  the  planes  associated  with  these  partitions 
are  denoted  by  0,,  i  —  l,...,n.  The  crease  edge  detection  algorithm  then  proceeds  as 
follows: 

(i)  All  points  with  G,  <  t,  t  l,...,n  are  discarded  for  some  t.  This  step  is  intended 
to  eliminate  deep  “interior”  points  that  are  surface  points  far  from  edges  or  jump 
boundaries. 

(ii)  The  likelihood  of  each  of  the  partitions  at  the  remaining  points  (denoted  by  set  11) 
is  computed  as  follows: 

II  7>(?/(*o,yo)|/i,),  (0 

(*<>  ,V<>  )*  < ! 
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where  p(tf(xo,  Vo)|-fy)  a  function  of  the  physical  properties  of  the  target  (surface 
orientation  with  respect  to  the  signal  beam  direction,  reflectance  and  actual  distance 
from  the  sensor)  and  sensor  parameters  (fixed  for  a  given  scene). 

Although  it  would  be  ideal  to  derive  p(ff(*o»yo)i^u)  from  the  physics  of  the  range  image 
function,  it  is  a  very  complicated  problem  that  is  not  currently  addressed.  Presently  we  are 
modelling  the  function  p(y(xo,yo))  as  normally  distributed  with  mean  n  =  distance  D  to 
the  estimated  plane  at  point  (xo,y<j)>  and  constant  variance.  If  reliable  reflectance  data  is 
available  then  the  variance  model  can  be  augmented  to  o(io*yo)  ~  <r(*0.*o)(^>P>  ^)>  where  p 
is  the  surface  reflectance  and  <f>  i  the  orientation  of  the  appropriate  plane  Rj  with  respect 
to  the  signed  beam  direction.  When  a  is  constant  the  likelihood  function  reduces  to  the 
root  mean  square  error  of  the  two  planes  combined. 

Once  the  likelihood  values  are  computed,  the  partition  with  maximum  likelihood  (or 
minimum  RMSE)  is  selected.  Then  all  points  for  which  fcK  <  t,  where  0  is  the  angle 
associated  with  the  best  partition,  are  identified  an  crease  edge  pouru 

A  more  implementation-oriented  description  of  the.  software  we  are  currently  using  for 
crease  edge  detection  is  given  here. 

1.  At  each  point  an  edge  with  a  particular  orientation  (horizontal  vertical,  either  di¬ 
agonal)  is  hypothesized.  Pixels  on  either  side  of  this  edge  fu.ti  to  two  separate 
planes. 

2.  The  root  mean  square  error  of  the  computed  planar  fit  I.  calculated. 

3.  When  the  fits  and  RMSEs  have  been  computed  for  each  of  *  e  edge  directions,  then 
the  edge  direction  with  minimum  RMSE  is  retained  as  a  potential  edge. 

4.  Potential  edges  are  then  accepted  or  rejected  based  upon  an  RMSE  threshold  and 
local  ainimum  filtering. 

5.  The  angle  between  the  two  surfaces  is  computed  for  each  of  the  accepted  edges.  If  this 
angle  is  less  than  a  certain  angular  threshold  then  the  edge  is  considered  nonexistent. 

This  technique  will  find  some  jump  edges  as  well  as  crease  edges  so  an  additional  test 
has  been  added  to  remove  jump  edges  based  upon  a  jump  distance  computed  between  the 
two  planes  at  the  point  in  question.  This  jump  distance  is  best  described  with  reference 
to  Figure  6.  The  point  where  the  jump  distance  Is  to  be  computed  is  lz.belled  x.  The 
planes  Py  and  P2  on  either  side  of  the  proposed  edge  cross  position  x  at  points  Zy  and  z2 
respectively.  The  jump  distance  is  the  average  of  the  shortest  distance  from  point  zj  to 
plane  P2  and  likewise  the  shortest  distance  from  z2  to  When  two  planes  are  parallel 
this  distance  is  identical  to  the  shortest  distance  between  the  planes.  This  measure  is  much 
more  dependent  on  relative  planar  orientation  than  on  absolute  orientation  which  makes 
it  more  useful  than  other  jump  distance  measures. 

An  example  of  the  behavior  of  crease  edge  defection  prior  to  thresholding  is  given  in 
Figure  7.  The  first  frame  is  the  base  image  that,  is  to  be  analyzed.  The  surface  segmen¬ 
tation  of  this  image  was  included  in  the  September  monthly  report.  Figure  7-b  shows  the 
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Figure  6:  Computing  the  Jump  Distance. 

minimum  RMSE  values  computed  at  each  pixel.  Bright  areaa  indicate  a  large  amount 
of  error  and  are  not  good  candidates  for  crease  edge  pixels.  Figure  7-c  shows  the  angles 
between  the  planes  of  best  fit.  Angles  which  are  too  small  indicate  no  crease  edge.  Finally 
the  last  frame  shows  the  jump  distance  corresponding  to  this  fit.  Large  jump  values  also 
eliminate  the  possibility  of  a  crease  edge.  Inspection  of  these  images  indicates  that  the 
most  likely  place  for  crease  edges  is  at  the  lower  edge  of  the  target,  between  the  body  of 
the  tank  and  the  turret,  and  within  the  turret  itself. 

Since  we  do  not  have  any  accurate  noise  models  for  the  laser  range  data  the  noise  model 
assumed  in  this  implementation  is  additive  Gaussian.  This  also  allows  for  easy  simplified 
testing  on  synthetic  data.  When  the  correct  noise  model  for  the  data  is  known,  it  can  be 
used  to  tailor  the  likelihood  function. 

When  crease  edge  detection  is  being  used  in  conjunction  with  the  surface  fitting  it  is 
not  necessary  to  test  for  crease  edges  over  the  entire  image.  The  crease  edge  detection 
can  be  restricted  to  those  areas  in  the  immediate  vicinity  of  surface  intersections  that 
result  from  the  surface  growing.  The  surface  growing  algorithms  also  use  thresholds  that 
are  related  to  the  likelihood  of  a  planar  fit.  These  thresholds  are  computed  initially  from 
those  portions  of  the  image  that  are  expected  to  contain  the  target.  The  threshold  used 
during  the  surface  growing  is  expected  to  be  highly  correlated  with,  if  not  identical  to,  the 
threshold  of  the  likelihood  function  mentioned  above. 
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Figure  7:  Crease  edge  example:  (a)  original  image,  (b)  error  for  the  plane  fit  associated 
with  best  edge  orientation,  (c)  angle  between  planes,  (d)  jump  distance  between  planes. 
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3.0  SUMMARY  AND  FUTURE  WORK 

Tha  progress  made  in  the  15  May  1987  to  15  November  1987  reporting  period  concerned 
preliminary  implementation  of  surface  segmentation/description  techniques.  A  modified 
version  of  Beal  arid  Jain’s  region  growing  segmentor  is  in  place  and  shows  promising  re¬ 
sults.  This  is  a  major  step  towards  developing  an  effective  LADAE  target  recognition 
methodology. 

In  the  previous  reporting  period  (November  1986  -  May  1987)  a  preliminary  simplified 
version  of  the  matching  algorithm  was  implemented.  Future  work  will  be  aimed  at  creating 
the  links  between  these  two  algorithms.  The  major  tasks  involved  in  this  effort  include: 

•  Integration  of  the  jump  and  crease  edge  detection  into  the  surfac  segmentation,  so 
that  a  final  segmentation  image  is  produced  composed  of  known  planes. 

•  Smoothing  and  multiscale  analysis  will  be  investigated  further  for  handling  the  ques¬ 
tion  of  variable  resolution. 

•  Description  stage:  Features  of  each  of  the  surfaces  and  their  relationships  must  be 
computed.  Any  sur  faces  thought  to  belong  to  the  background  must  be  removed  at 
this  stage.  The  feasibility  and  utility  of  curve  description  will  be  investigated  as  well. 

•  Matching:  The  matching  routines  must  be  upgraded  to  handle  a  larger  feature  set, 
occlusion,  and  varying  ranges.  Orientation  approximation  routines  will  be  added  to 
eliminate  the  need  for  multiple  target  views  in  the  database.  The  utility  of  different 
features  as  an  aid  to  matching  will  be  studied. 

•  Database:  A  search  must  be  made  for  imagery  that  more  accurately  approximates 
the  expected  working  data.  A  major  problem  in  testing  has  been  the  scarcity  of 
data  representative  of  the  specifications  that  were  originally  assumed.  It  would  be 
useful  to  be  able  to  run  tests  on  data  which  is  free  of  range  ambiguities  and  is  not 
range-gated. 

•  Modelling:  Graph  models  for  several  target  types  will  be  developed  based  upon  the 
ERIM  wireframe  models. 
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